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Abstract. We introduce and analyze two-level and multi-level preconditioners for a family 
of Interior Penalty (IP) discontinuous Galerkin (DG) discretizations of second order elliptic 
I problems with large jumps in the diffusion coefficient. Our approach to IPDG-type methods is 

based on a splitting of the DG space into two components that are orthogonal in the energy 
^SJ inner product naturally induced by the methods. As a result, the methods and their analysis 

— depend in a crucial way on the diffusion coefficient of the problem. The analysis of the pro- 
^ posed preconditioners is presented for both symmetric and non-symmetric IP schemes; dealing 

simultaneously with the jump in the diffusion coefficient and the non-nested character of the 
relevant discrete spaces presents extra difficulties in the analysis which precludes a simple ex- 
tension of existing results. However, we are able to establish robustness (with respect to the 
diffusion coefficient) and nearly-optimality (up to a logarithmic term depending on the mesh 
size) for both two-level and BPX-type preconditioners. Following the analysis, we present a 
sequence of detailed numerical results which verify the theory and illustrate the performance of 
the methods. The paper includes an Appendix with a collection of proofs of several technical 
results required for the analysis. 



1. Introduction 

^ Let C IR'^ be a bounded polygon (for d = 2) or polyhedron (for c? = 3) and / G L^(0). We 

consider the following second order elliptic equation with strongly discontinuous coefficients: 

CN| / -V-(kVu) = / inO, 

\ u = Q ondn. ^ ' 

The scalar function k = k{x) denotes the diffusion coefficient which is assumed to be piecewise 
constant with respect to an initial non-overlapping (open) subdomain partition of the domain 
^ n, denoted Ts = {Qm}m=i^ with U^^^Tlm = H and n r2„ = for n / m. Although the 

^ (polygonal or polyhedral) regions ,m = 1 . . . M, might have complicated geometry, we will 

always assume that there is an initial shape-regular triangulation To such that = k{x)\t 



is a constant for all T G To- Problem (1.1) belongs to the class of interface or transmission 



X 

u 

problems, which are relevant to many applications such as groundwater flow, electromagnetics 
and semiconductor device modeling. The coefficients in these applications might have large 
discontinuities across the interfaces between different regions with different material properties. 



Finite element discretizations of (1.1) lead to linear systems with badly conditioned stiffness 
matrices. The condition numbers of these matrices depend not only on the mesh size, but also 
on the largest jump in the coefficients. 

Much research has been devoted to developing efficient and robust preconditioners for con- 



forming finite element discretizations of (1.1). Nonoverlapping domain decomposition precondi- 
tioners, such as Balancing Neumann-Neumann [48j, FETI-DP [45] and Bramble-Pasciak-Schatz 
Preconditioners [12] have been shown to be robust with respect to coefficient variations and 
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mesh size (up to a logarithmic factor), in theory and in practice, but only if special exotic 
coarse solvers (such as those based on discrete harmonic extensions fSTl HH [38] ) are used (see 
also [60] ) . The construction and use of such exotic coarse spaces is avoided in other multilevel 
methods, such as the Bramble-Pasciak-Xu (BPX) or multigrid preconditioners, for which it has 
always been observed that when used with conjugate gradient (CG) iteration, result in robust 
and efficient algorithms with respect to jumps in the coefficients, independently of the problem 
dimension. However, their analysis (based on the standard CG theory) predict a deterioration 
in the rate of convergence with respect to both the coefficients and the mesh size. By resorting 
to more sophisticated CG theory (see [Gj Section 13.2], [7J) which accounts for and exploits 
the particular spectral structure of the preconditioned system^ the authors in |58l [61] show 
that standard multilevel and overlapping domain decomposition methods lead to nearly optimal 
preconditioners for CG algorithms. (See also [26]). Much less attention has been devoted to non- 
conforming approximations. Overlapping preconditioners for the lowest order Crouzeix-Raviart 
approximation of (1.1) are found in [52, 51j, where the analysis depends on the assumption that 
the coefficient k is quasi-monotone. 

In this article, we consider the construction and analysis of preconditioners for the Interior 
Penalty (IP) Discontinuous Galerkin (DC) approximation of ( 1.1 ). Based on discontinuous finite 
element spaces, DG methods can deal robustly with partial differential equations of almost any 
kind, as well as with equations whose type changes within the computational domain. They are 
naturally suited for multi-physics applications, and for problems with highly varying material 
properties, such as ( 1.1 ). The design of efficient solvers for DG discretizations has been pursued 
only in the last ten years; and, while classical approaches have been successfully extended to 
second order elliptic problems, the discontinuous nature of the underlying finite element spaces 
has motivated the creation of new techniques to develop solvers. Additive Schwarz methods 
(of overlapping and non-overlapping type) are considered and analyzed in [39l IMl [21 131 [U [H] . 
Multigrid methods are studied in [HI [201 [El [13 IMl 129]. Two-level methods are presented 
in [311 '22', "23] . More general multi- level methods based on algebraic techniques are considered 
in |47l i46j . However, all the analysis in these works consider only the case of a smoothly or slowly 
varying diffusivity coefficient. For problem ( 1.1 ), only in |34l 1351 136] have the authors introduced 



and analyzed non-overlapping BBDC and FETI-DP domain decomposition preconditioners for 
a Nitsche type method where a Symmetric Interior Penalty DG discretization is used (only) on 
the skeleton of the subdomain partition, while a standard conforming approximation is used 
in the interior of the subdomains. Robustness and quasi-optimality is shown in d = 2 for the 
Additive and Hybrid BBDC [35J and FETI-DP [36j preconditioners, even for the case of non- 
matching grids. As it happens for conforming discretizations, the construction and analysis of 
these preconditioners rely on the use of exotic coarse solvers, which might complicate the actual 
implementation of the method. 

The goal of this article is to design, and provide a rigorous analysis of, a simple multilevel 
solver for the lowest order (i.e. piecewise linear discontinuous) approximation of a family of 
Interior Penalty (IPDG) methods. To ease the presentation, we focus on a minor variant of the 
classical IP methods, penalizing only the mean value of the jumps: the "weakly penalized" or 
IPDG-0 methods (called Type-0 in [lOj). Our approach follows the ideas in ^Oj, and it is based 
on a splitting of the DG space into two components that are orthogonal in the energy inner 
product naturally induced by the IPDG-0 methods. 

Roughly speaking, the construction amounts to identifying a "low frequency" space (the 
Crouzeix-Raviart elements) and then defining a complementary space. However, a notable 
difference takes place in the DG space decomposition introduced for the Laplace equation 



Namely, that there are a few small eigenvalues due to the jump coefficient distribution that have no influence 
in the (observed) overall convergence of the iteration 
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[To] . For problem (1.1), the subspaces depend on the coefficient k, and this is certainly related 



to the splittings used in algebraic multigrid (AMG jl5j). With the orthogonal splitting of 



the DG space at hand, the solution of problem (1.1) reduces to solving two sub-problems: a 



non-conforming approximation to (1.1), and a problem in the complementary space containing 
high oscillatory error components. We show the latter subproblem is easy to solve, since it is 
spectrally equivalent to its diagonal form, and so CG with a diagonal preconditioner is a uniform 
and robust solver. 

For the former subproblem, following \58\ I61j. we develop and analyze (in the standard and 
asymptotic convergence regimes) a two-level method and a BPX preconditioner. Nevertheless, 
dealing simultaneously with the jump in the coefficient k and the non-nested character of the 
Crouziex-Raviart (CR) spaces presents extra difficulties in the analysis which precludes a simple 
extension of |58| I61j . We are able to establish nearly optimal convergence and robustness 
(with respect to both the mesh size and the coefficient k) for the two-level method and for 
the BPX preconditioner (up to a logarithmic term depending on the mesh size). The resulting 
algorithms involve the use of a solver in the CR space that is reduced to a smoothing step 
followed by a conforming solver. Therefore, in particular one can argue that any of the robust 



and efficient solvers designed for conforming approximations of problem ( 1.1 ) could be used as a 
preconditioner here. Finally we mention that, although the two- level and multilevel methods we 
propose are based on the piecewise linear IP-0 methods, they could be used as preconditioners 
for the solution of the linear systems arising from high order DG methods. 

Outline of the paper. The rest of the paper is organized as follows. We introduce the 
IPDG-1 and IPDG-0 methods for approximating (1.1) in ^and revise some of their properties. 



The space decomposition of DG finite element space is introduced in q3} Consequences of the 
space splitting are described in Q The two-level and multi-level methods for the Crouzeix- 
Raviart approximation of (1.1) are constructed and analyzed in ^ Numerical experiments are 



included in f|6j to verify the theory and assess the performance and robustness of the proposed 
preconditioners. In ^ we briefly comment on how the developed solvers and theory can be 
extended for the classical IPDG-1 family. The paper is completed with an Appendix where we 
have collected proofs of several technical results required in our analysis. 

Throughout the paper we shall use the standard notation for Sobolev spaces and their norms. 
We will use the notation xi < yi, and X2 > y2i whenever there exist constants Ci, C2 indepen- 
dent of the mesh size h and the coefficient k or other parameters that xi, X2, yi and y2 may 
depend on, and such that xi < Ciyi and X2 > 6*22/2 j respectively. We also use the notation 
j; ~ y for Cix < y < C2X. 



2. Discontinuous Galerkin Methods 
In this section, we introduce the basic notation and describe the DG methods we consider 



for approximating the problem (1.1) 



Let Th he a shape-regular family of partitions of O into d-simplices T (triangles in d = 2 
or tetrahedra in d = 3). We denote by Ht the diameter of T and we set h = maxT^Th ^T- 
We also assume that the decomposition Th is conforming in the sense that it does not contain 
hanging nodes and that Th C To, with To being quasi- uniform initial triangulation that resolves 
the coefficient k. We denote by £h the set of all edges/faces and by and the collection of 
all interior and boundary edges/faces, respectively. The space H^{Th) is the set of element-wise 
functions, and L'^{£h) refers to the set of functions whose traces on the elements of £h are 
square integrable. 

Following [5], we recall the usual DG analysis tools. Let and T~ be two neighboring 
elements, and let n"*", n~ be their outward normal unit vectors, respectively (n^ = ny±). Let 
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and be the restriction of C and r to T^. We set: 

2|C}} = (C+ + r ), ICl = C+n+ + C n- on e G 4°, 

2{{r}} = (r+ + r-), [r] = t+ • n+ + ■ n" on e G S"^. 

We also define the weighted average, §-§5, for any 6 = {6e}ee£° with G [0, 1] Ve: 

= 5eC+ + (1 - <5e)C- , iT}s,=6eT+ + {l-6e)T- , On e G 4° . (2.1) 



For e G f^, we set 

[Cl=Cn, |r}} = {{T}k=T oneG4^. (2.2) 
We will also use the notation 

{u,'w)'ji= / uwdx M u,w £ L^{Q), {u,w)£^= / uwds \/ u,w,£ L'^{£h)- 

The DG approximation to the model problem can be written as 

Find ^ G V^^ such that ^^^(n^^, tu) = (/, w)T^, , Vu; G V^^ , 

where V^*^ is the piecewise linear discontinuous finite element space, and A^'-'{-, •) is the bilinear 
form defining the method. 

In this paper, we focus on a family of weighted Interior Penalty methods (see |54|). with special 
attention given to a variant (weakly penalized) of them. The bilinear form defining the classical 
family of weighted IP methods [54|, here called IP(/3)-l methods, is given by A^'^{-, •) = A{-, •), 
with 

A{v,w) = {KyhV,Vw)T^ - {iKVv'^p^,lw\)e^+e{[vliKVw}p:)£^ 

+ {aK^Kelvl H)f„ V^;, w G V^"^ . ^^'^^ 

where = -I gives the SIPG(/3)-l methods; 9 = 1 leads to NIPG(/3)-l methods; and 6* = 
gives the IIPG(/3)-l methods. Here, he denotes the {d — 1) dimensional Lebesgue measure of 
e G f/i.The penalty parameter a > is set to be a positive constant; and it has to be taken large 
enough to ensure coercivity of the corresponding bilinear forms when 9^1. The symmetric 
method was first considered in [5l] and later in [33^ Section 4] for jump coefficient problems 
(although there it was written using a slightly different notation and DG was only used in the 
skeleton of the partition) . It was later extended to advection-diffusion problems in ^25j and [30j . 

We also introduce the corresponding family of IP(/3)-0 methods, which use the mid-point 
quadrature rule for computing the integrals in the last term in (2.3) above. That is, we set 
^^^(•,-)=A(-,-) with 

+ {ah-^KeV%v\), H)£„ Vr;, G V^"" , 

where '■ L'^i^h) ^ ^^i^h) is the L^-projection onto the piecewise constants on S^. We 
note that this projection satisfies H'Pe = 1. In (2.3) and (2.4), for any e G f ^ with 

e = 5T+ n dT~ , the coefficient kt and the weight fie are defined as follows: 

k.t = k\t, I3e = ^_ , where = k\^^, (2.5) 

The coefficient the harmonic mean of and k : 

Ke := ^ _ . 2.6 

The weight /3 = {/3e}ee£'° depends on the coefficient k and therefore it might vary over all 
interior edges/faces (of the subdomain partition To resolving the coefficient k). 



MULTILEVEL PRECONDITIONERS FOR DG METHODS FOR JUMP COEFFICIENT PROBLEMS 



Remark 2.1. We note that one could take Ke as mm{K,~^ , k }, since both are equivalent: 



min {k , K } < Kg 



— < 2min{K^,K } < 2k 



± 



(2.7) 



The equivalence relations in (2.7) show that the results on spectral equivalence and uniform 



preconditioning given later for (2.3) with Kg defined in ( |2.6[ ) (the harmonic mean) will auto- 
matically hold for method (2.3) with := mm{K~^ , To fix the notation and simplify the 
presentation, we stick to definition (2.6) for Kg. 

Weighted Residual Formulation. Following |21) we can rewrite the two families of IP meth- 
ods in the weighted residual framework: For all v,w £ V^*^, 

A{v,w) = (-V • {KVv),w)r, + {l^'^vliw^i_p^)so + {lvlBi{w))e^, 

A,{v,w) = (-V • {KVv),w)r, + {{^^yvlMi-p.)s^, + (H,^°(^iH)>£., 

where Bi is defined as: 



(2.8) 
(2.9) 

(2.10) 



Throughout the paper both the weighted residual formulation (2.8)-(2.9) and the standard one 
(2.3)-(2.4) will be used interchangeably. 

We now establish a result that guarantees the spectral equivalence between A{-^ ■) and Ao{-, •). 



Lemma 2.2. Let A{-,-) he a bilinear form corresponding to a IP(/3)-l method (2.3) and let 
Ao{-,-) be the corresponding IP((3)-0 bilinear form as defined in (2.4). Then there exists a 
positive constant cq = co(a), depending only on the shape regularity of the mesh and the penalty 
parameter a (but independent of the coefficient k and the mesh size h ) such that, 



Ao{v, v) < A{v, v) < co{a)Ao{v, v) £ 



DG 



(2.11) 



Proof. The lower bound follows immediately from the fact that the projection is an L'^{Sh)- 
orthogonal projection and therefore has unit norm. The upper bound would follow if we show 



^ ah~' KeWHWl, <C{Y, ^tIMIt + ^K^^eWVl 



which can be proved by arguing exactly as in [101 [13 E] and taking into account (2.7). 



□ 



By virtue of Lemma |2.2[ it will be enough throughout the rest of the paper to focus on 
the design and analysis of multilevel preconditioners for the IP(/3)-0 methods. At least in the 
symmetric case, the preconditioners proposed for SIPG(/3)-0 will exhibit the same convergence 
(asymptotically) when applied to SIPG(/3)-l. 



Continuity and Coercivity of IP(/3)-0 methods. The family of methods ( 2.4 ) can be shown 
to provide an accurate and robust approximation to the solution of (1.1). We define the energy 
norm: 



\v\ 



DGO ■— / , 



^ ^t||Vz;||2 + V Kg/ig-i||PgO(H)||2 



(2.12) 



Then, ^o('; ■) is continuous and coercive in the above norm, with constants independent of the 
mesh size h and the coefficient k: 



Continuity: 
Coercivity: 



\Aq{v,w)\ < \\\v\ 
•^o{v,v) > \\\v\ 



DGO 



w 



I DGO ' 



(2.13) 
(2.14) 
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Although the proof of (2.14) and (2.13) is standard, we sketch it here for completeness. Note 



first that for each e £ £^ such that e = dT^ n dT , the weighted average -{[kVi^J^^ can be 
rewritten as: 

U^v}}^^ = /3e(K+(Vr;)+) + (1 - (3e){^^- (Vv)-) 



K 



+ K~ 



(2.15) 



Trace inequality [Tj, inverse inequality [28j and (2.7) imply the following bounds 



< 2(Ke)a(l + CD (k+\\Vv\\It++K-\\Vv\\It- 



This inequality, combined with Cauchy-Schwarz inequality and (2.7), gives 



< 



-c J e 
-t-h 

L. — ^ rv 



1/2 



1/2 



Now (2.13) follows from Cauchy-Schwarz inequality. The inequality (2.14) is proved by setting 
w = V in ( |2.3"1 ) and taking into account the above estimate. We have then 

Aoiv,v) = E '^t||V^;||2_^ + a E ^^e/ie-'ll^e (H)llo,e " (1 " ^)({{^Vt;}U, H)^, 



Ten 



(H)ll 



0,e ' 



and (2.14) follows immediately by taking a > 1 large enough (if 9 1). Moreover, notice that 
both constants in (2.13) and (2.14) depend on the shape regularity of the mesh partition but 



are independent of the coefficient k. 



Obviously, continuity and coercivity also hold for the IP(/3)-l methods (2.3) if the norm (2.12) 
is replaced by 



\dg ■= E '^2^ll^^llo,r + E ^^^e ^1 



IO,e- 



(2.16) 



See |33j or [9] for a detailed proof. For both families of methods, optimal error estimates in 
the energy norms (2.12) and (2.16) can be shown, arguing as in [5]. See also f8] for further 
discussion on the L^-error analysis of these methods. 



MULTILEVEL PRECONDITIONERS FOR DG METHODS FOR JUMP COEFFICIENT PROBLEMS 



7 



3. Space decomposition of the V,^'^ space 
In this section, we introduce a decomposition of the V^^'^-space that wih play a key role in 



the design of the solvers for the DG discretizations (2.3) and (2.4). In \W\ 124) . it is shown 
that the discontinuous piecewise linear finite element space V^*^ admits the decomposition: 
V^'^ = ® Z, where denotes the standard Crouzeix-Raviart space defined as 

V^^ ={vG L\n) : vy G F^T) \fT G % and • n) = Ve G , (3.1) 

and the complementary space Z is a space of piece-wise linear functions with average zero at 
the mass centers of the internal edges/faces: 

Z = {z£ L^{Q) : G F^{T) VT G % and 7^°(l^l) = 0, Ve G f^} . 

In [lU], it was shown that this decomposition satisfies Aq{v^ z) = when k = 1, for all v G 
and z & Z. We now modify the definition of Z above in order to account for the presence of a 



coefficient in the problem (1.1). Let 



Zp = {ze L\n) : z\^ G ¥\T) VT G % and P°({{^>i-/3e) = 0, Ve G £t] , (3.2) 



where the weight /3e was defined earlier in (2.5). Note that the weight /3e depends on the coef- 



ficient K, and, as a consequence, the space Z/^ is also coefficient dependent. In what follows, we 
shall show that Zjs is a space complementary to in V^*^ and the corresponding decompo- 
sition has properties analogous to the properties of the decomposition = © Z given 
in [10] for the Poisson problem. 

For any e £ £h with e C T G 7/t, let ipe^T be the canonical Crouzeix-Raviart basis function on 
T, which is defined by 

<^e,T|T e P^(T), ipe^Tim^') = 4,6' Vc' G £h(T), and ife^rix) = Vx T, 

where rue is the mass center of e. We will denote by ut and ue the number of simplices and 
faces (or edges when d = 2) respectively. We also denote by n^E the number of boundary faces. 

Proposition 3.1. For any u G V^*^ there exists a unique v G V^'^ and a unique zp G Zp such 
that u = V + zp , that is 

y^"" = V^''®Z^. (3.3) 

Proof. For simplicity, throughout the proof we will set /?"'" = /3e, /3~ = (1 — /3e), and ip^ = <fe,T± 
for any e G with e = dT^ n dT~ . We also denote = fe,T for any e G with e = dTDdi}. 
Since the mesh is made of d-simplices 

dim Vj^'^ = {d+ l)nT = 2nE — nsE, 
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and it is also obvious that {^t}eG£° U We}e£e^ form a basis for Vj^'-' . Notice that /3+ + = 1, 
we can therefore express any u G as 

u{x) = ^ u^{me 

eesi ee£° ee^e 

= iP'^^i^e) + /3+n-(me))(^+(x) + <f-{x)) 

+ ^{U^i^e) - u"(me))(/3+V3+(3;) - P~'ip~{x)) + ^ u(me)93e(x) 



= w(x) + Z/3(x). 

Then for each e G f^, we set 

:=(/p+(x) + <^-(x), (3.4) 



and ipe{x) := for all x T+UT . In the definition (3.5) of il)l{x), we have used ^p^ (x) = for 
X G T"*" and V5^(x) = for x G . Finally, when e G with e = 5T n dVL for some T, we set 

ijl(x) = ^e(x), VxGT. (3.6) 
It is then straightforward to check that 

= span{v?f^}ee£-°, and Z/j = spanlV'ejee^-ft- 
Hence, for all u G Vy^^ there exist unique v G Vy^^ and G defined by 

rCR 



such that u = V + zp. This shows (3.3) and concludes the proof. □ 

Remark 3.2. As we pointed out in the introduction, the definition of the subspace clearly 
depends on the coefficient k, since /3 depends on k. Such dependence is often also seen in 
algebraic multigrid analysis, where the coarse spaces depend on the operator at hand. They are 
in fact explicitly constructed in this way, the aim being to increase robustness of the methods. 



In the proof of Proposition 



3.1 



above, we have introduced the basis in both V^^ and Zj^. 



The canonical Crouzeix-Raviart basis functions jeefo are continuous at the mass centers 
me of the faces e G 8'^. The basis {•i/'fjeefh ™ consists of piecewise functions, which are 
discontinuous across the faces in 8^. In fact, for any z ^ Zp such that z = X^eef^ •^eV'l with 
G IR, we have 

(|zjn+)(me') = Ze', Ve' G 8h. 
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To see this, evaluating the jump of z at gives 

(Hn+)(m,0 = ^Ze(|Ve1n+)(meO = v([V'e'In+)(meO 

^e'(/?e'-(/3e'-l))=^e', e' E f^, 

This relation will also be used later to obtain uniform diagonal preconditioners for the restric- 
tions of ^(•, •) and .4o(-, •) on Z^. 

Remark 3.3. For mixed boundary value problems, that is, dil. contains both Neumann boundary 
Tat 7^ and Dirichlet boundary T o with dVt = T/j U F^r , the definition of the basis functions on 



the boundary faces [see (3.6)] needs to be changed as: 



</,f«(x) = ^Pe,T{x), e = dTnTN, for all x E T, , . 

^l(x) = (peA^i e = dTnTD, for all x E T. ^ ' 

Thus, in case Ftv 7^ the dimension of is increased (by adding to it functions that 

correspond to degrees of freedom on F^r) and the dimension of is decreased accordingly. 
Clearly things balance out correctly: the identity = V^^Q) Z/s holds, and also the analysis 
carries over with very little modification. 

Next lemma is a simple but key observation used in the design of efficient solvers. 



Lemma 3.4. LetAo{-,-) be the bilinear form defined in (2.4). Then 



Ao{v,z) = VwEVf^, VzEZ^. (3. 



Furthermore ifAo{-,-) is symmetric (and positive definite) then the decomposition (3.3) is Aq- 
orthogonal, namely, _L^(, Z/^. 



Proof. From the weighted-residual form of Ao{-, •) given in for ah V E V^^, and all z € Z^ 



we easily obtain 

Ao{v,z) = (-V • {KVv),z)r, + (I^V^l,{{^li-/3j^;j + {lvlV^emz)))e, = 0. 

In the equation above, the first term is zero due to the fact that v E V^^j so v is linear in 
each T, and the coefficient k E P'^(r). Last term vanishes (independently of the choice of 9, 
or equivalently the choice of 13i{v)), because {M,V^{Bi{z)))£^ = 0, thanks to the definition 
(3.1) of V^^. The second term vanishes from the definition of Zjj (since I^Vf] is constant on 



each e E <S^). Moreover, in the case when ^o(') ") is symmetric and positive definite we have 
that Ao{v,z) = Ao{z,v), for all v E and for all z E Zjs. Thus, for the symmetric method 
^o(") ■)) the spaces and Zj^ are indeed ^o-orthogonal. The proof is complete. □ 

4. Solvers for IP(/3)-0 methods 



In this section we show how Proposition |3.1| and Lemma |3.4| can be used in the design and 
analysis of uniformly convergent iterative methods for the IP(/3)-0 methods. We follow the ideas 
and analysis introduced in [10] and point out the differences. We first consider the approximation 
to problem ( |1.1[ ) with A^'^{-, •) = Ao{-, •). To begin, let ^0 be the discrete operator defined by 



{Aou,w) = Ao{u,w) and let Aq be its matrix representation in the new basis (3.4) and (3.5). 
We denote by u = [z, v]-^, f = [fz, fv]"^ the vector representation of the unknown function u and 
of the right hand side /, respectively, in this new basis. A simple consequence of Lemma 3.4 is 
that the matrix Aq (in this basis) has block lower triangular structure: 

Aq^ 
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where Aq^jA™ are the matrix representation of Aq restricted to the subspaces Zjj and V^^, 
respectively, and Ag^ is the matrix representation of the term that accounts for the couphng 
(or non-symmetry) Ao{tl^^, ip^^). As remarked earher, for SIPG(/3)-0, the stiffness matrix Aq is 
block-diagonal. 

gives a 2D example, with two squares = [—0.5,0]^ and U2 = [0,0.5]^ inside the 
[—1, 1]^. We set the coefficients hi{x) = 1 for all x £ QiU Q2 and k{x) = 10~^ for 



4.1 



Figure 
domain Cl 

X G Q, \ (Oi U $12)- Figures 4.2 and 4.3 show the sparsity patterns of the IP(/3)-0 methods with 





Figure 4.1. Computational domain and unstructured mesh. 



standard nodal basis and the basis (3.4)-(3.5), respectively. 



SIPG (nnz- 10834} 




NIPG (nnz=10834) 
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Figure 4.2. Non-zero pattern of the matrix representation in the standard 
nodal basis of the operators associated with IP(/3)-0 methods. From left to 
right: SIPG, NIPG and IIPG methods. 

Clearly, as in the constant coefficient simple algorithm based on a block version 

of forward substitution provides an exact solver for the solution of the linear systems with 
coefficient matrix Aq. A formal description of this block forward substitution is given as the 
next Algorithm. 

Algorithm 4.1 (Block Forward Substitution). 

1. Find z £ such that Aoizj'ip) = {f,i^)Th all ip £ Zjs 

2. Find v £ V^^ such that Ao{v, tp) = (/, ip)r,, - Ao{z, ip) for all p G V^^ 

3. Set u = z + V 

The above algorithm requires the solution of ^o(") ") on Z^ (Step 1. of the algorithm) and 
the solution of ^o(-, •) on (Step 2. of the algorithm). Unlike the situation in ^Oj, due to 



MULTILEVEL PRECONDITIONERS FOR DG METHODS FOR JUMP COEFFICIENT PROBLEMS 11 




Figure 4.3. Non-zero pattern of the matrix representations according to the 
basis spHtting (3.3) of the operator associated with IP(/3)-0 methods, i.e., Aq. 
From left to right: SIPG, NIPG and IIPG methods. 

the jump coefficient in (1.1), the solution on is more involved, and therefore we postpone 
its discussion and analysis until Section [5j We next discuss the solution on Zj^. 

4.1. Solution on Zp. In this section we describe the main properties of the IP(/3)-0 methods 
when restricted to the Zp, which will in turn indicate how the solution of Step 1. of Algorithm 



4.1 can be efficiently done. 



The first result in this subsection establishes the symmetry of the restrictions of the bilinear 
forms of the IP(/3)-0 methods to Zp. 



Lemma 4.2. Let Aq{-,-) be the bilinear form of a IP(fi)-0 method as defined in (2.4). Then, 
the restriction to Z/^ of Ao{-, ■) is symmetric. Namely, for 9 = —1, 0, 1, we have 

Proof. If 6 = —1 there is nothing to prove, since in this case both bilinear forms are symmetric. 
Hence we only consider the cases 6 = or 9 = —1. Integrating by parts and using the fact that 
z G Zi3 and G Zp are linear on each element T shows that 

= (-V • (kVV^), Vz)r, = (^V^, Vz)r, - (|/^VV^}}/3., lz})£, - {lnVnAz}}i-p.)eo . 

Hence, from the definition (3.2) of the Zj^ space, it follows that 

iKVi;,Vz)r, = iUVij}},,^, lzj)e, = {ifiVz}}^^, m)e„ Vz, V e Z^. (4.2) 

Substituting the above identity in the definition of the bilinear form (2.4) then leads to 

Ao{z,i;) = 9{lzliKVn^.)s, + (^°(W),'^eM)^, 

= 9{KV4^,Vz)r, + {V°em),^^elz})£,=Ao{%b,z). 

This shows the symmetry of ^o("i ") on Zjj, regardless the value of 9. □ 

We now study the conditioning of the bilinear form Ao{-, •) on Zjj. For all z G Z/^, and for all 
(f) G Zi3 with 

Z = ^ Zeipe ^ ^13, and = ^ (l^elpe ^ ^P- 

we introduce a weighted scalar product (•, •)=„ : Zjj x Zjj ^ M, and the corresponding norm || • ||*, 
defined as follows 

{Z,(p)^ := "^^-^KeZeCpe : \\z\\l := [z, z)^. (4.3) 



Observe that the matrix representation (in the basis given in ( |3.5[ ) ) of the above weighted scalar 
product is in fact a diagonal matrix. The next result shows that the restriction of Ao{-, ■) to Zjj 
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is spectrally equivalent to the weighted scalar product (•,•)* defined in (4.3) and therefore its 
matrix representation Aq^ is spectrally equivalent to a diagonal matrix. 



(4.4) 



Lemma 4.3. Let he the space defined in ( |3.2[ ). Then, the following estimates hold 

Ml<Mz,z)<\\z\\l -izeZp. 



Proof. Let us fix z E -Z^, z = YleeSh "^e^l- From the definition of "Ped-^D' is immediate to 
see that 



l^e°(W)llo,e 



\e\z1. 



Thus, we have that 



^IDIlL 



E 



Zip. 



(4.5) 



To show the upper bound in (4.4), we notice that (4.2) together with (2.7) and the standard 
trace and inverse inequalities gives 

E ^t\\Vz\\It = {KVz,Vz)r, = {{{kVz}^^, lzj)e^ = (^e|Vz|,pO(H))£, 



TeTh 



1/2 



1/2 



< I ht\\Vz\\1^t 



zlDIlL 



and therefore by (14. 5|), 



IO,e 



(4.6) 



Since z ^ Zfs was arbitrary, we have that Ao{z, z) < ||z||^ for all z G Zf^. This proves the upper 
bound in (4.4). 

To prove the lower bound, we use the coercivity estimate (2.14) for the bilinear form ^o('; ■) 
in the energy norm ||H||£)co [^^^ ( 2.12[ )]. For all z G Zjs we have 



Aoiz,z) 



> 



z 



IO,e 



\z\ 



which is the desired bound and gives (4.4). 



□ 



Last result guarantees that the linear systems on Zj^ can be efficiently solved by preconditioned 
CG (PCG) with a diagonal preconditioner. As a corollary of the result in Lemma 4.3, the 
number of PCG iterations will be independent of both the mesh size and the variations in the 
PDE coefficient. 

We end this section by showing that in the particular case of the IIPG(/3)-0 method, the 
matrix representation of ^o(") ") restricted to is in fact a diagonal matrix. (See the rightmost 
figure in Fig. 4.3). 



Lemma 4.4. Let .4o(-,-) be the bilinear form of the non-symmetric LIPG{f3)-0 method (2.4) 
with 9 = 0. Let {V'f }ee£^h be the basis for the space Z^ as defined in (3.5). Let Ag^ be the matrix 
representation in this basis of the restriction to the subspace Zp of the operator associated to 
Ai){-,-). Then, Aq^ is diagonal. 
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Proof. Note that from the definition ( |2.4[ ) of the method {6 = 0) together with (4.2) we have 



(4.7) 



Let {■(/'|}ee£-^ be the basis functions (3.5). To prove that Aq^ is diagonal it is enough to show 
that for the basis functions (3.5), the following relation holds: 



(4. 



where Se,e' is the delta function associated with the edge/face e. We now show (4.8). Observe 
that the supports of and ip^, have empty intersection unless e, e' C T for some T £ Th- Let 
Tn dQ = be an interior element, then from (|4.7|) and the mid-point integration rule, we have 



ah-' / KeV^,im)V^eiire4)ds = ah-'Ke[2re{me)]m'{rne)] 
J e 

Aah~' Ke6e,e' , 6,6 CdT, 6,6 G£^, 



which shows (4.8) for interior edges with Ce 



4a/i„ 



For boundary edges/faces the con- 



siderations are essentially the same and therefore omitted. The proof is complete since the 
relation (4.8) readily implies that the off-diagonal terms of Ag^ are zero. □ 



5. Robust Preconditioner on V^^ 
In this section, we develop efficient and robust (additive) two-level and multilevel precondi- 



tioners for the solution of the IP(/3)-0 methods in the CR space (cf. Step 2 of algorithm 4.1 ). We 



first review a few preliminaries and tools that will be needed for the convergence analysis. We 
then define the two- level preconditioner and provide the convergence analysis. The last part of 
the section contains the construction and convergence analysis of the multilevel preconditioner. 

From the definition (3.1) of the space, it follows that the restriction of ^o(') ■) to 
reduces to the classical P^-nonconforming finite element discretization of (1.1): 

FindnGVf^: Ao{u,w) = {KVu,Vw)r,, = if,w), WweV^^. (5.1) 

We denote as the operator induced by ( |5.1[ ). For the analysis in this section, we will need 
the following semi- norms and norms for any v G Vh^'- 



|2 

l0,K 



M 



\0,T ' 



\i,h,n, ■- 



E 



\Vv\ 



0,T 5 



l0,K 



+ 



\'^\l,h,K 



(5.2) 



(5.3) 



Since (5.1) is a symmetric and coercive problem, from the classical theory of PCG we know 
that the convergence rates of the iterative method for Aq^ with preconditioner, say B, are 
fully determined, in the worst case scenario, by the condition number of the preconditioned 
system: }C{BAq^). However, if the spectrum of BAq'^, a{BAQ^) happens to be divided in 
two sets: a{BA^^) = ao{BA^^) U ai{BA^^), where ao{BA^^) = {Ai, . . . , A^} contains all of 
the very small (often referred to as "bad") eigenvalues, and the remaining eigenvalues (bounded 
above and below) are contained in ai{BAQ^) = {Am+i, ■ • • , A^^^}, that is, \j G [a,b] for 
j = m + 1, . . . , ncR, with ncR = dim(y^'^) = ue — nBE, i-e. the number of interior edges, then 
the error at the fc-th iteration of the PCG algorithm is bounded by (see e.g. [6t Hij 17]): 



k—m 



\U ■ 



Uk\\l,h,K 
■ 'Uo\\i^h,K 



< 2{IC{BA^ 



CR\ 



(5.4) 
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The above estimate indicates that if m is not large (there are only a few very small eigenvalues) 
then the asymptotic convergence rate of the resulting PCG method will be dominated by the 



factor ^'^=Jt-^, i.e. by ^/b/a where b = X]\[{BAq^) and a = Xm+i{BAQ^). The quantity (b/a) 

which determines the asymptotic convergence rate is often called effective condition number. 
This is precisely the situation in the case of problems with large jumps in the coefficient k. 



In fact, for a conforming FE approximation to (1.1) it has been observed in |42^ [58] that the 
spectrum a{BAQ^) might contain a few very small eigenvalues, which result in an extremely 
large value of )C{BAq^). Nevertheless, they seem to have very little influence on the efficiency 
and overall performance of the PCG method. Therefore, it is natural to study the asymptotic 
convergence in this case, which as mentioned above is determined by the effective condition 
number: 

Definition 5.1. Let be a real A/^-dimensional Hilbert space, and ^ :V V he a symmetric 
positive definite linear operator, with eigenvalues < Ai < • • • < Aat. The m-th effective 
condition number of 21 is defined by 

Ajv(2l) 



A„^+l(2l)• 



Below, we will introduce the two-level and multilevel preconditioners, and study in detail 
the spectrum of the preconditioned systems. In particular, we gave estimates on both condition 
numbers and the effective condition numbers, which indicates the pre- asymptotic and asymptotic 



convergence rates in (5.4) of the PCG algorithms 



5.1. Two-level preconditioner for ^o(") ") on V^^. In this subsection, we construct a two- 
level additive preconditioner, which consists of a standard pointwise smoother ( Jacobi, or Gauss- 
Seidel) on the nonconforming space plus a coarse solver on a (possibly coarser) conforming 
space V~°^^ '■= {v € ^^o(^^) : v\t G IPiC^), VT € 7^}. Here, T~ refers to a possibly coarser 

partition such that 7^ ^ 7^; that is for h = h, Tj^ is the same as Th, while for h = H > h 
the partitions are nested and Th could be regarded as a refinement of 7r. Observe that V~°^^ 

h 

is a proper subspace of Vf^^. To define the two-level preconditioner, we consider the following 
(overlapping) space decomposition of Vf^^: 

^ = + . (5.5) 



On V-'^'^^ we consider the standard conforming P^-approximation to (1.1): Find x ^ V~"^^ 



h 

such that 



h 



Aoix,v)=aix,r]) = J^f^^X-V7jdx={f,r^), Vr? G . (5.6) 



The bilinear form in (5.6) defines a natural "energy" inner product, and induces the following 
weighted energy norm: 

\x\l.,D-= [ ^\^x\^dx, yx^H^D), Den. (5.7) 
Jd 

For simplicity, we write |x|i,k = |x|i,K,r2 denote by A^' the operator associated to (5.6). We 



define the two level preconditioner as: 

B : ^ ^ ^, B := R'^ + {A^y^Q^, (5.8) 

where is the operator corresponding to a Jacobi or symmetric Gauss-Seidel smoother on 
V^^, and Q'^ : i— t- is the standard L-^-projection. We refer to [9] for further details 

on the matrix representation of the above preconditioner. 



MULTILEVEL PRECONDITIONERS FOR DG METHODS FOR JUMP COEFFICIENT PROBLEMS 15 



Next Theorem is the main result of this section, which estabhshes the convergence for the 



two- level preconditioner (5.8). 



Theorem 5.2. Let B be the multilevel preconditioner defined in ( |5.8[ ), and w = h/h he the 

ratio of the mesh sizes of andTh- Then, the condition number K,{BAq^) satisfies: 



(5.9) 



where := maxx^Th ^t/ ^^^TeTh is what we refer as the jump of the coefficient and 

Co > is a constant independent of the coefficient k and the mesh size. Moreover, there exists 
an integer rriQ depending only on the distribution of the coefficient k such that the mo-th effective 
condition number lCmo{BAQ^) satisfies: 



where Ci > is a constant independent of the coefficient and mesh size. Hence, the convergence 
rate of the PCG algorithm can be bounded as 



< 



\u 



Uo\l,h,K 



2{CoJiK)vj^log{2w)-iy 



rCxw\og^l'^{2vj) - \ 
fCxvj\oi'l'^(2w) + 1 



k—mg 



(5.10) 



Remark 5.3. We emphasize that for two-level preconditioners, since the ratio w = h/h is a fixed 
constant, the effective condition number }Cmo{BAQ^) is bounded uniformly with respect to the 
coefficient variation and mesh size. Clearly, according to estimate (5.10), the number of (pre- 
asymptotic) PCG iterations will depend on the constant mo (the number of floating subdomains; 
see (5.19)). While such a bound could be a large overestimate (depending on the coefficient 
distribution), it is sufficient for our purposes. Since rriQ is fixed, the asymptotic convergence 
rate in (5.10) is bounded uniformly with respect to coefficient variation and mesh size. In short, 



while the estimates given here might not be sharp with regard to the pre-asymptotic PCG 
convergence, they are asymptotically uniform with respect to the parameters of interest. 



We recall the following well known identity [59, Lemma 2.4]: 

X,v-x) + aix,x)] V?;GVf^, 



{B-^v,v)= inf [n{v 



(5.11) 



where TZ^-,-) is the bilinear form associated with the smoother defined by TZ{v,w) := {Rv,w) 
for any uj,v ^ ^h '^- "^^^ proof of Theorem 5.2 amounts to showing a smoothing property for 
7^(-, ■) and the stability of the decomposition given in (5.5 ). Next Lemma establishes the former; 
the latter is contained in next subsection. 

Lemma 5.4. LetTl{-,-) be the bilinear form associated to Jacobi, or symmetric Gauss-Seidel 
smoother. Then we have the following estimates 



•^q{v,v) <TZ{v,v) and TZ{v,v)'^h 



ll^ll0,K ) 



(5.12) 



Proof. We only need to show this inequality for Jacobi smoother, since the Jacobi and the sym- 
metric Gauss-Seidel methods are equivalent for any SPD matrix, see for example |56|, Proposi- 
tion 6.12] or |631 Lemma 3.3]. 

For any v E V'^f^, we write v = See^^ '^eH^'t^ where ip^^ is the basis function with respect 
to e G Note that for Jacobi smoother, we have 



CR ^CR^ 
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For any e G let £{e) := {e' e £^ : e' C dT, T e Th dT D e}. Then, Cauchy-Schwarz 
and the arithmetic-geometric inequahties give 

Ao{v,v) = Yl MV>e^,V>e'^)VeVe' 

ee£-° e'e£{e) 

= E M^e"", ^eVe < E My^e"", ^eVe = CsTliv, v). 

e&£° e^£l 

The constant Cs above only depends on the cardinality j^E[e)^ which is bounded by 5 in 2D and 
7 in 3D. This proves the first inequality in (5.12). 



Since the mesh is quasi-uniform, for any v = Vef^ ^ ^^'^ ^ have 



eCdT 



Now by direct calculation, for any basis function (f^^ we have 



0,K,T 



CRii2 



\0,K,T- 



Therefore, by the equivalence relations (5.13) and ( |5.14 ), we get 



l0,K 



E ^'livv^: 



CR\\2 

0,K,T+UT- 



(5.13) 



(5.14) 



E E -eiiv^f«iit.T- E E h~'<y^'wl.,T 

TeTh eCdT TeTh eCdT 

TeTh 



which concludes the proof. 



□ 



5.2. A stable Decomposition. In this subsection we give a detailed discussion of the sta- 
ble decomposition. The main tool is an operator : — t- ^"'^^ that satisfies certain 
approximation and stability properties, as stated in the next Lemma. 



Lemma 5.5. There exists an interpolation operator P^ : 



that satisfies the fol- 



CR 



lowing approximation and stability properties: 

Approximation: || (/ - P^)i;||o,k < Cah\log2h/h\^/'^\\v\\i^h,K, 

Stability: < log2h/h\'/^\\v\\i,h,. G V^^, 

with constants Ca and Cs independent of the coefficient k and mesh size. 

A construction of such an operator P^, and proof of the above results, are given in the Appen 



(5.15) 
(5.16) 



dix 



A 



We would like to point out that the operator P/* is not used in the actual implementation 



of the preconditioner B, as it is plainly seen from (5.8). However, the operator and its 
approximation and stability properties play a crucial role in the analysis. 

Observe that on the right hand side of (5.15) and (5.16), the bounds are given in terms of the 
weighted full iJ^-norm ||f In general, one cannot replace the norm by the energy 
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norm induced by the bilinear form ah{-, •)• To replace the full norm by the semi-norm, 

one might use the Poincare-Friedrichs inequality for the nonconforming finite element space (cf. 
[321 US]) to get: 



V 



,2 ^ f \ [ \ \2a <: ( ^ I |2 ^ maxreT^^KT, ,2 

lo,K < / 1^1 dx< max^r |t;|i^fe< l^liA^- 



From the above inequality, we have: 

........ ;^ _ ^ 



Corollary 5.6. There exists an interpolation operator P!^ : Vf'^ — )■ V-"'^^ satisfying the follow- 



ing approximation and stability properties: 

\\{I-Phv\\o,. < j'/\K)h\log{2h/h)\'/^\vh,H,. , yveV^"", (5.17) 
\Ph\i,>^ < j'/\K)\log{2h/h)\'/MiA- , yveV^'', (5.18) 
where J{k) = max^^gT^ ht/ minTeTh '^r jump of the coefficient. 



The approximation and stability properties given in Corollary 5.6 depend on the coefficient 
variation However, by imposing some constraints on the finite element space V^^, it is 

possible to get rid of this dependence obtaining a robust result. Following |55l Definition 4.1] 
we introduce the index set of floating subdomains (the subdomains not touching the Dirichlet 
boundary) : 

T:={i : measd_i(5Sl n 5X7^) = } . (5.19) 
We then introduce the subspace C V^^: 

V^R ■= \ v£ ^ : / vdx = Wiel]. (5.20) 

The key feature of the above subspace is the fact that the Poincare-Friedrichs inequality for 
nonconforming finite elements space (cf. \32\ 116) ) now holds on each subdomain, which allows 
us to replace the full norm ||f by the semi-norm |f for any v G V^^- 

We remark that the condition on the zero-average in (5.20), is not essential; other conditions 
could be used (see [55|) as long as they allow for the application of a Poincare-type inequality. 
At this point, we would like to emphasize that the dimension of V^'^ is related to the number 
of floating subdomains and in fact: dim(V^^) = dim(V'^^) — rriQ, where tuq = #X is the 
cardinality of I. 

By restricting now the action of the operator Pj^ to functions in V^^, we have the following 
result, as an easy corollary from Lemma 5.5 Its proof follows (as mentioned above) by applying 
Poincare-Friederichs inequality (for nonconforming) on each subdomain. 

Corollary 5.7. Let C he the subspace defined in (5.20). Then, there exist an operator 
Ph ■ ^h^ ^f^^ satisfying 

||(/-Pj)i-||o,« < h\\og{2h/h)\'/Mi,h,>^. yv^V^^, (5.21) 

\pIv\i,k < \log{2h/h)\'/^\v\,,h,. , Vt;GVf«. (5.22) 

With the aid of the results from Corollary |5.6| and Corollary |5.7[ we can finally show the 
stability of the decomposition (5.5). 

Lemma 5.8. j 

property holds: 



Lemma 5.8. For any v G Vh^' ^ ~ Phi'") ^ ^^^^ following stable decomposition 



n{v - x,v - x) + aix,x) < Ji^^)ih/hy\log2h/h\\v\l,^^^ . (5.23) 
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In particular, for any v G we have 



n{v -x,v-x) + a(x, X) < {h/hf\ log 2h/h\ 



(5.24) 



Proof. Below, we give a proof ( 5.24[ ). Given any v G V^^, let x £ V~ be defined as x := Ph^. 



By the approximation property (5.21) of given in Corollary 



5.7 



we have 



Xllo,re 



h''^\\v 



PML<{h/hr\log2h/h\\v\i,,, 



where in the first inequality, we have used (5.12) from Lemma 5.4 For the second term, the 



stability (5.22) of Pt' from Corollary 



5.7 



gives, 



aix,x) = \Piv\l,<\log2h/h\\v\lf^^^. 
The proof of (5.24) is complete. The proof of ( 5.23| ) is essentially the same but using Corol- 



lary instead of Corollary 5.7 



□ 



We have now all ingredients to complete the proof of Theorem 5.2 



Proof of Theorem 5.2. To estimate the maximum eigenvalue of BAq^, let x ^ V~^^^ and v G 

be arbitrary. We set vq = {v — x), and so v = vq + x- The Cauchy-Schwarz inequality and 
Lemma |5.4| yield 

Aq{v, v) = Aq{vq + x,vq + x) < 2(A(i'o, vq) + Aq{x, x)) < d {T^{vo, vq) + a{x, x)) , 



where ci = 2max{cs, 1}, with Cg (defined in the proof of Lemma 5.4), is a constant independent 
of K and mesh size. Using the identity (5.11) and the fact that x £ arbitrary, we have 



Hence, 



^max(^^o'^) 



Ao{v,v) <ci{B-^v,v), \fveV^^. 
= max — — j = max —— r < ci. 



which is uniformly bounded, independently of the coefficient and the mesh size. 

Let w = h/h be the ratio of the mesh sizes. For the lower bounds of Amin and Amo+ii 
with X = Ph'^ together with (5.11) give 



Lemma 



{B-^v, v) <n{v-x,v-x) + \X\1 . < J{t^)^^\ log 2w\Ao{v, v), Vt; G ^, 



[B ^v, v) < Tl{v — x^v 
The first inequality implies that 



x) + IxiL < ^'1 \og2w\Ao{v,v), yv G ^. 



mm 



Aoiv,v) 



> 



vev^R {B-^v,v) ~ J{K)w'^\log2w\' 



The second inequality, together with the fact that dim(V'^^) = dim(V^'^) — mo and the 
minimax principle [40, Theorem 8.1.2]) gives 



A™o+i(S<^) > niin 



An{v,v) 



> 



gy-Cfl (B ^V,v) ^ w'^\\og2-!u\ 

Therefore, the condition number )C{BAq^) and the effective condition JCmoiBA^^) can be 
respectively bounded by 

1C{BA^'^) < CoJ{k)tu^\ log2CT|, and ICmoiBA^^) < Cxw^\ \og2w\, 

with Co and Ci, constants independent of the coefficient and mesh size. The inequality (5.10) 
then follows directly from (5.4). □ 
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5.3. Multilevel Preconditioner for Ao{-,-) on . We now introduce a multilevel pre- 



conditioner, using the two-level theory developed before. The idea is to replace [A'-''] ^ in (5.8) 
with a spectrally equivalent operator B*-^ : V-""^^ i— )• V~°^^, which corresponds to the additive 
BPX preconditioner (see e.g. [HI [57]). 

Given a sequence of quasi-uniform triangulations T^- for j = 0, 1, • • • , J, we denote by Wj = 
y^onf (^j — Q^i^ . . . ^ J") g^Tij(^ consider the family of nested conforming spaces (defined w.r.t. the 

family of partitions {T^I^^q): 

WoCWiC ■■■ CWj . 

Here, we assume that the coarsest triangulation To resolves the jump in the coefficient, and 
without loss of generality, we also assume that hj ~ (j = 0, • • • , J) and h = hj. The space 
decomposition that we use to define the multilevel BPX preconditioner is: 



CR 



j=0 



J+1 

j=0 



(5.25) 



where we have denoted VFj+i = V^'^. For j = 0, • • • , J we denote by A'j^ the operator corre- 
sponding to the restriction of a(-, •) to Wj, namely 

{A^Vj,Wj) = a{vj,Wj), \/vj E Wj, \/wj S Wj. 
The operator form of the multilevel preconditioner then reads: 



J+i 



Bml : Vf ^ ^ 



Bml ■■= [A^r^Q^ + Rj'Qj- (5-26) 

Here, Qj : i— )• Wj is the L2-orthogonal projection on Wj for j = 0, . . . , J and we set 
Qj+i = I- We use an exact solver on the coarsest grid. With this notation in hand, one can 
prove that 

J+i 

a{wo,wo) + 'ETZj{wj,Wj) . (5-27) 



ML^' 



V,V 



E 



inf 



1, . . . , (J + 1) correspond to Jacobi or symmetric Gauss-Seidel smoothers, and 



Here Tlji-, •), 

the proof of (5.27) is similar to (5.11) for the two-level case. 



Next two results will be used in our convergence analysis. 

Lemma 5.9 ( |58[ Lemma 4.2]). Let lZj{-,-) he the Jacobi or the symmetric Gauss-Siedel 
smoother for the solution of the discretization (5.6) on Wj space (ij = 1, • • • , J). Then, 



a{w,w) < nj{w,w) < h- ^||w||o,«; Vw G Wj. 
We also need the following strengthened Cauchy Schwarz inequality. 

Lemma 5.10 (Strengthened Cauchy Schwarz, cf. j57| Lemma 6.2]). For j = 1, ■ 

j < I < J , there exists a constant 7 G (0, 1) such that 

a{wi,Wj) < 7'"-'(V^II^HIo,K)(^7^ll^illo,K), ^wi G Wi,Wj G Wj. 
The main result of this section is the following: 



,J — 1 and 



(5.28) 



Theorem 5.11. Let Bml be the multilevel preconditioner defined in (5.26). Then, the condition 
number )C{BmlAq'^) satisfies: 

/C(Sml<'^) < CoJ(k)J2 , 
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where J is the number of levels, and := maxj^gT^ kt/ minj^gT^ kt is the jump of the 

coefficient. Moreover, there exists an integer mg depending only on the distribution of the 
coefficient k such that the m^-th effective condition number }Cmo{BML^o^) satisfies: 

where the constants Co,Ci > are independent of the coefficients and mesh size. Hence, the 
convergence rate of the PCG algorithm can be bounded as 



< 2{CqJ{k)j^ - 1; 



mo 



(5.29) 



Proof. We first give a bound on Amax(-BML^o'^)- Let v S VFj+i = Vf^^ be arbitrary, and let 
{wj}jt,o be any decomposition of v, namely v = X]/=o ""^j; with Wj G Wj. By the Cauchy- 
Schwarz inequality, we have 



J+i J+i 



J J 



Ao{v,v) 



i=i j=i 



By Lemma [5. 9[ the strengthened Cauchy-Schwarz inequality (Lemma 5.10), and the smoothing 
property of 7^J+l(•, •) ( |5.12[ ), we get: 

J J 



Ao{v,v) 



< 



a{wo,wo^ 



+ '''' ^11 



=1 i=i 
J 



< a{wo,wo) + ^TZj{wj,Wj) +TZ{wj+i,wj+i] 
\ j=i 

where in the second inequality, we used the fact that the spectral radius of the matrix {■y^^~^^)jxj 
is uniformly bounded by (1 — 7)"^. Since the decomposition of v was arbitrary, taking the 
infimum above over all such decompositions and using the identity (5.27) then gives 

Aoiv,v)<iB^lv,v), V^;Gy;f^, 

which shows that Amax(-SML^o'^) ^ 1- 



Similar to the proof of Theorem 5.2, the estimates on the lower bound for Amm and A. 



mo 



rely on the stability of the decomposition. For this purpose, we make use of the interpolation 
operator and its properties introduced in ^5.2 To simplify the notation, we set Pj := P^^ : 

V^^ — ^ for j = 0, . . . , J, and set Pj+i = I and P_i = 0. Given any v G V^^, we define 
the decomposition of v as 

J+i J+i 
V = Pj+iv = ^^(-Pj — Pj-i)v = Wj, where vuj = {Pj — Pj-i)v. 

j=0 j=0 

Clearly, Wj G Wj for j = 1, • • • , ( J + 1) and wq = Pqv £ Wq. Triangle inequality and the 
smoothing properties of TZj {j = 1, - ■ ■ , J + 1) from Lemma 5.9 and Lemma 5.4, give 

J+i J+i 
aiwo,wo) + Y.nj{wj,Wj) < \Pov\l, + hfUPj " Pj-iHln 
j=i i=i 

J 

j=0 
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Using in (5.30), the approximation property (5.17) of Pj (j = 0, • • • , J) and the stabihty property 



(5.18) of Pq given in Corohary 5.6, we obtain 

J+i 

(-^ML^'^) ^ a{wo,wo) + 

i=i 

This gives the estimate on minimal eigenvalue of -Bml^o"^ as 



CR 



Similarly, if we use Corollary 5.7 in (5.30), we obtain 



{B^lv,v)<j'Ao{v,v), Vr;eyf^. 
Therefore, Xmo+iiB^AQ^) > 1/J^ by the minimax principle and the result follows. 



□ 



Remark 5.12. Similar results hold also for the multiplicative multilevel methods such as the 
y-cycle. These results can be derived from estimates comparing multiplicative and additive 
preconditioners given in [431 Theorem 4] or |27l Theorem 4.2]. We refer to [62] for a detailed 
analysis and numerical justification. 

6. Numerical Experiments 



We consider the model problem (1.1) in the square Q, = [—1, 1]^ with coefficients: 

-0.5,0]2 U [0,0.5]2, 



k{x) 



1.0, 
e, 



Vx G [- 
elsewhere. 



In all of the following experiments, e varies from 10~^ up to 10^, covering a wide range of 
variations of the coefficients. The set of experiments is carried out on a family of structured 
triangulations; we consider uniform refinement with a structured initial triangulation on level 
with 32 elements and mesh size h = 2^^. This initial mesh resolves the jump in the coefficients. 
Each refined triangulation is then obtained by subdividing each element of the previous level 
into four congruent elements. The number of degrees of freedom Ni in the DG discretizations 
on each level satisfies Ni = A^Nq for £ = 0, 1,2,3,4 with A^o = 96. We consider the IP(/3)-0 



method (2.4) with penalty parameter a 



We use the basis ( |3.4[ )-( |3.5[ ) for the computations 
use Algorithm 14.1 



To solve the resulting linear systems we 
Due to the block structure(4.1 ) of Aq (matrix representation of Aq in the 



basis (3.4)-(3.5)) we only need to numerically verify the effectiveness of the solvers for each 

is the same (since 
is 



block; A^^ and Ag^ 



Recall that for any choice of 9 = 0, ±1, 



the block A™ 



it is the stiffness matrix of the Crouzeix-Raviart discretization (5.1)), while the block Aq 
different for different values of 9, but it is always an SPD matrix. To solve each of these smaller 
systems we use the preconditioned CG, for which we have set the tolerance to TOL=10~^ for the 
stopping criteria based on the residual; namely, if tq is the initial residual and r^ is the residual 
at iteration k, the PCG iteration process is terminated at iteration k if ||?^'^||£2/lk'^lk2 < 10"'^. 
The experiments were carried out on an IMAC (OS X) with 2.93 GHz Intel Core 17, and 8 GB 
1333 MHz DDR3. 

The systems corresponding to Ag^ are solved by a PCG algorithm using its diagonal as 
a preconditioner. The estimated condition numbers of D~^Aq^ for SIPG(/3)-0 are reported in 
Table 6.1 Observe that the condition numbers of 
1, which confirms the result established in Lemma 4.3 
its diagonal. Similar results, although not reported here, were found for the NIPG(/3)-0 and 



-1 A 22 
2 '^O 

J^Ag^ are uniformly bounded and close to 
i.e., that Aq^ is spectrally equivalent to 
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e 


levels 


h 




10-^ 


10-^ 


1 


10^ 


10^ 


lO'' 







1.73 (14) 


1.73 (12) 


1.73 


1.73 (9) 


1.72 (10) 


1.73 (12) 


1.73 (13) 


1 




1.72 (15) 


1.72 (13) 


1.72 


1.72 (10) 


1.72 (10) 


1.72 (12) 


1.72 (14) 


2 


2-' 


1.72 (15) 


1.72 (13) 


1.72 


1.71 (10) 


1.7 (10) 


1.71 (12) 


1.72 (15) 


3 


2-4 


1.72 (15) 


1.72 (12) 


1.71 


1.71 (10) 


1.69 (10) 


1.69 (12) 


1.69 (16) 



Table 6.1. Estimated condition numbers /C(D^ ""^Aq^) (number of PCG itera- 
tions) for the block Aq^ in SIPG(/3)-0 discretization. 



IIPG(/3)-0 methods. The system Aq'' arising from the restriction of ^o(")") to the Crouzeix- 



Raviart space is solved by a PCG algorithm with the two-level preconditioners defined in (5.8), 
for which we use 5 symmetric Gauss-Seidel iterations as smoother. 

Figure 6.1 shows the spectrum of the preconditioned system for e = 10~^ and the mesh size 
!~^. In this example, we have taken h = h, so Ij^ = Th- Note that there is only one (very 



h 




Figure 6.1. Eigenvalue distribution of 



I'' for e 



10-5 and h = 2- 



small) eigenvalue close to zero (which may be related to the fact that there are only 2 different 
values for the coefficients). In Table 6.2 we report the estimated condition number /C(BA™) and 
the effective condition number (denoted by /Ci(BAg'')). Observe that the estimated condition 
number /C(BAq^) deteriorates with respect to the magnitude of the jump in the coefficient. In 
contrast, the effective condition number /Ci(BAq") is uniformly bounded with respect to both 



the mesh size and the jump of coefficient, as predicted by Theorem 5.2 



For comparison, we also pres ent t he r esults obtained with different choices of coarse grid 

As we can see from these two tables, the effective 



6.3 



6.4 



h = 2h, Ah, reported in Tables 
condition number is uniformly bounded with respect to the coefficient and mesh size. However, 



comparing to the results in Table 6.2, it seems that the effective condition numbers get larger 



when we use a coarser grid. These observations coincide with the conclusion in Theorem 5.2 



(5.26). 



We now present the results corresponding to the multilevel preconditioners as defined in 
In Table 6.5 we report the estimated condition number /C(BAq^) and the effective 
condition number (denoted by /Ci(BAq'')) for the BPX. Also for the BPX, we use 5 symmetric 
Gauss-Siedel iterations as a smoother. Observe that the estimated condition number /C(BAq^) 
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e 


levels 





1 


2 


3 


4 


h 




2-^ 


2-3 


2-4 




IQ-^ 


/Ci(BAg^) 


4.52 


Q '^1 p_|_4 n q^i 
3.37 


9 77p-I-4 (I'i) 
2.95 


2 S7p-I-4 ('91 
2.78 


9 n8p-l-4 ('91 
2.71 


10-3 


IC (MAli'") 
/Ci(]BAg^) 


301 (11) 
4.48 


333 ns") 

3.36 


280 ('18'! 
2.95 


240 fl8l 
2.77 


211 (18) 
2.71 


10-1 


/Ci(]BA^^) 


4 49 n fTl 
2.97 


5 99 n 3^1 

fj . I XO J 

2.89 


4 91 (^ A\ 

'i. t? X 1 X4: 1 

2.69 


4 7 n 4^1 

t: . 1 \ X'X J 

2.6 


4 59 (^^\ 

t:. U t? I X'i J 

2.57 


1 


/Ci(]BAg^) 


2.16 (8) 
2.06 


2 25 fill 
2.16 


2 29 ('121 
2.21 


2.3 (12) 
2.19 


2 33 ('12') 
2.18 


10^ 


/C(]BAg^) 
/Ci(]BAg^) 


2.33 (9) 
2.3 


3.16 (12) 
2.63 


3.58 (13) 
2.66 


3.8 (14) 
2.62 


3.95 (14) 
2.61 


10^ 


/C(BAg^) 
/Ci(BA™) 


2.54 (9) 
2.4 


4.12 (13) 
2.82 


5.37 (14) 
2.85 


6.56 (15) 
2.8 


7.79 (16) 
2.78 


10^ 


/C(BAg^) 
/Ci(BAg'') 


2.55 (9) 
2.4 


4.13 (13) 
2.83 


5.41 (15) 
2.85 


6.62 (16) 
2.8 


7.89 (17) 
2.78 



Table 6.2. Two level preconditioner for Aq'' on with h = h. 



e 


levels 





1 


2 


3 


4 




h 


2-1 


2-' 


2-' 


2-4 


2-^ 


10-5 


/C(BAg^) 


X 


4.92e+4 (18) 


4.28e+4 (24) 


3.66e+4 (26) 


3.21e+4 (27) 




X 


4.27 


3.61 


3.38 


3.33 


10-3 


/C(BA^^) 
/Ci(BAg-") 


X 
X 


494 (16) 
4.26 


431 (21) 
3.61 


370 (21) 
3.38 


325 (21) 
3.34 


10-1 


/C(BA™) 
/Ci(BAJ^^) 


X 
X 


7.14 (14) 
3.46 


6.69 (16) 
3.27 


6.35 (16) 
3.2 


6.19 (16) 
3.19 


1 


/C(BA^^) 
/Ci(BAg^) 


X 
X 


2.63 (11) 
2.32 


2.75 (13) 
2.61 


2.91 (14) 
2.63 


2.97 (14) 
2.61 


lOi 


/C(BA^^) 
/Ci(BAJ^^) 


X 
X 


3.74 (13) 
3.33 


4.3 (15) 
3.38 


4.48 (16) 
3.32 


4.67 (16) 
3.29 


103 


/C(BA^^) 
/Ci(BA^^) 


X 
X 


4.93 (14) 
3.64 


6.59 (16) 
3.65 


8.02 (18) 
3.56 


9.55 (18) 
3.49 


10^ 


/C(BAS^) 
/Ci(BAJ^'') 


X 
X 


4.95 (14) 
3.65 


6.63 (16) 
3.65 


8.02 (18) 
3.53 


9.66 (19) 
3.49 



Table 6.3. Two level preconditioner for A™ on with h = 2h. 



deteriorates with respect to the magnitude of the jump in coefficient. On the other hand, the 
effective condition number /Ci(BA™) is nearly uniformly bounded with respect to both the mesh 
size and the jump of the coefficient, as predicted by Theorem 5.11 Moreover, we also observe 



that the effective condition numbers grow linearly with respect to the number of levels, which 
is better than the quadratic growth in Theorem 5.11 This issue will be further investigated in 
the future. 



7. Solvers for IP(/3)-1 Methods 

We now briefly discuss how the preconditioners developed here for the IP(/3)-0 can be used 
or extended for preconditioning the IP(/3)-l methods (2.3). We follow ^lOj . 
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e 


levels 





1 


2 


3 


4 


h 


2-1 


2-^ 


2-3 


2-4 


2-^ 


w-^ 




X 
X 


X 
X 


7 RQp-i-A ('SI 

6.58 


7 9Qp-i-4 
5.99 


fi 41 p_|_4 Cit^) 
5.97 


10-3 


ICi{MA^^) 


X 
X 


X 
X 


6.57 


733 (28^1 
5.99 


646 f29l 
5.97 


10-1 


/Ci(]BAJj^) 


X 
X 


X 
X 


1 9 9 ('9n'l 
5.58 


1 1 6 l'99'l 
5.69 


1 1 4 ('99') 
5.76 


1 


IC (MA}i'") 
/Ci(]BAJj^) 


X 
X 


X 
X 


4.73 (17) 
3.99 


5 22 fl9'l 
4.75 


5 32 ('19') 
4.8 


10^ 


/C(]BA^^) 
/Ci(]BA;J'") 


X 
X 


X 
X 


7.55 (19) 
6.34 


6.84 (21) 
5.63 


6.97 (22) 
5.95 


10^ 


/C(]BA^^) 
/Ci(]BA;J^) 


X 
X 


X 
X 


11.2 (20) 
6.99 


12.2 (23) 
6.11 


14.6 (25) 
6.39 


10^ 


/C(BAg^) 
/Ci(BA;j^) 


X 
X 


X 
X 


11.3 (20) 
7 


12.3 (23) 
6.12 


14.9 (26) 
6.4 



Table 6.4. Two level preconditioner for Aq'' on V^^ with h = 4/i. 



e 


levels 





1 


2 


3 


4 


h 




2-^ 


2-3 


2-4 


2-^ 


10-5 


/C(BA^'') 
/Ci(BAS'') 


3c+4 (12) 
4.52 


5.03e+4 (27) 
5.69 


6.77e+4 (33) 
6.81 


8.64e+4 (37) 
7.9 


1.06C+5 (42) 
9.03 


10-3 


/C(BAS") 
/Ci(BA™) 


301 (11) 
4.49 


506 (22) 
5.65 


680 (27) 
6.78 


868 (31) 
7.86 


1.06e+03 (35) 
8.98 


10-1 


/C(BAJJ^) 
A:i(BAS'') 


4.42 (10) 
2.97 


7.5 (16) 
4.22 


9.92 (20) 
5.28 


12.5 (24) 
6.3 


15.1 (26) 
7.41 


1 


/C(BA^'') 

a:i(bas^) 


2.16 (8) 
2.07 


3.32 (13) 
3.17 


4.45 (17) 
4.25 


5.61 (20) 
5.23 


6.67 (22) 
6.24 


lOi 


/C(BAS^) 

a:i(ba™) 


2.33 (9) 
2.3 


4.58 (14) 
3.84 


6.69 (19) 
5.06 


8.75 (22) 
6.19 


11 (26) 
7.31 


103 


/C(BAS^) 
/Ci(BAS^) 


2.54 (9) 
2.4 


5.92 (16) 
4.11 


10.1 (21) 
5.42 


15.6 (25) 
6.62 


23 (29) 
7.81 


10^ 


/C(BAS'') 
/Ci(BAS'') 


2.55 (9) 
2.4 


5.94 (16) 
4.11 


10.2 (21) 
5.43 


15.7 (25) 
6.62 


23.3 (29) 
7.81 



Table 6.5. PCG with BPX (additive) preconditioner for solving on 



7.1. Solvers for the SIPG(/3)-l method. From the spectral equivalence given in Lemma 



2.2, it follows that any of the preconditioners designed for ^o("; ■) result in an efficient solver for 



A{-, •). Motivated by the block diagonal form of Aq (cf. ( |4.1| )), we use the decomposition (3.3) 
and define the following block-Jacobi preconditioner: 

Block-Jacobi: B^^ := [R^]-^ + BQ^^ , (7.1) 

where denotes the operator corresponding to the diagonal of A{-,-) restricted to and 
B refers to the corresponding multilevel preconditioner for the symmetric SIPG(/3)-l method 



(i.e., including the jump-jump term). The next result is a simple consequence of Theorem 5.11 



(focusing only on the asymptotic result) together with Lemma 2.2 
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Theorem 7.1. Let be the preconditioner defined through (7.1). Let niQ be the number of 

floating subdomains. Then, the following estimate holds for the effective condition )Cmo (B^^A) : 



The constant C > above is independent of the variation in the coefficients and mesh size. 



In Table 



7.1 



are given the estimated condition numbers of /C(]Bf"^A) together with the es- 
timated effective condition numbers /Ci(]Bf"^A), and the number of PCG iterations required 
for convergence. As can be seen from these two tables, /C(Bf"^A) deteriorate rapidly when e 



e 


levels 





1 


2 


3 


h 




2-^ 


2-3 


2-^ 


IQ-^ 


/Ci(BfGA) 


2.85C+4 (44) 
6.27 


3.37e+4 (44) 
6.33 


3.1e+4 (46) 
6.45 


2.85e+4 (46) 
6.49 




/C(Bf'^A) 
/Ci(BfGA) 


288 (33) 
6.24 


340 (34) 
6.3 


313 (34) 
6.42 


289 (32) 
6.46 


10-1 


/C(Bf'^A) 
/Ci(Bf^A) 


7.25 (22) 
5.62 


7.33 (22) 
5.6 


7.21 (22) 
5.71 


7.13 (22) 
5.73 


1 


/C(Bf'^A) 
/Ci(BfGA) 


5.53 (19) 
5.17 


5.76 (20) 
5.45 


5.8 (20) 
5.46 


5.83 (20) 
5.46 


lOi 


/C(Bf'^A) 
/Ci(BfGA) 


6.66 (22) 
5.91 


7.16 (23) 
6.2 


7.16 (23) 
6.25 


7.43 (23) 
6.27 


10^ 


/C(Bf'^A) 
/Ci(BfGA) 


6.38 (27) 
5.51 


8.98 (30) 
6.53 


11.1 (31) 
6.59 


13.5 (32) 
6.59 


10^ 


A:(Bf'^A) 
/Ci(BfGA) 


6.91 (33) 
6.38 


9.02 (36) 
6.54 


11.3 (39) 
6.6 


13.8 (40) 
6.59 



Table 7.1. Estimated condition number /C(Bf"^A) (number of PCG iterations) 



and the effective condition number /Ci(Bf A). 



becomes smaller, but /Ci(Bf"^A) are nearly uniformly bounded with respect to the coefficients 



and mesh size. These results confirm the theory predicted by Theorem 7.1 



7.2. Solvers for the non-symmetric IIPG(/3)-l and NIPG(/3)-l methods. We consider 
the following linear iteration: 

Algorithm 7.2. Given initial guess uq, for /c = 0, 1 . . . until convergence: 

1. Set ek = BDG{f -Auk); 

2. Update n^+i = Uk + e^ ■ 

Here, A : i— )• is the operator associated with the bilinear form of either the 

NIPG(/3)-l or IIPG(/3)-l methods with 6* = 1 and (9 = 0, respectively): 

{Av, w) := A{v, w), Vv, w G V^^ . (7.2) 

Following [lOj we consider as preconditioner B^'^ the symmetric part of A, defined by: 

B^^:=A^\ where {Asv,w):=^[A{v,w)+A{w,v)], \fv e V^^ , Vu; € F/f ^. (7.3) 

We note that from this definition and ( 2.14[ ), we immediately have that As is symmetric and 
positive d efini te. The next result guarantees uniform convergence of the linear iteration in 
Algorithm 7.2 with preconditioner B^'^ given by (7.3). The proof follows [10\ Theorem 5.1] and 
it is omitted. 
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Theorem 7.3. Let a* be a fixed value of the penalty parameter for which the IIPG(/3 )-0 bilinear 
form (2.4) ^o(') ") is coercive. LetA{-,-) be the bilinear form of the IIPG(fi)-l method ( 2.3|) 
with penalty parameter a > 4a*. Let B^^ = A^^ be the iterator in the linear iteration 7.2, 
and let Uk and ut+i be two consecutive iterates obtained via this algorithm. Then there exists a 
positive constant A < 1 such that 



\u ■ 



■"fc+llllDG ^ ^ 111^ ~ '<J-k\\\DG 



(7.4) 



To verify Theorem 
to the IIHII^x^) 



7.3 



we have computed the ^-norm (which is obviously equivalent in Vj^^ 
of the error propagation operator: E = I — B^^A = I — A'Z^A, for different 



meshes and values of e. This norm gives us the contraction number of the linear iteration in 



Algorithm 7.2, and so an estimate for the constant A in Theorem 7.3 The results are reported 
in Table El 







e 


levels 


h 


10"^ 


10-^ 


10"^ 


10"^ 


10-^ 


1 


10^ 


10^ 


10^ 


10^ 


10^ 







0.20 


0.20 


0.20 


0.20 


0.20 


0.20 


0.19 


0.19 


0.19 


0.19 


0.19 


1 


2-^ 


0.14 


0.14 


0.14 


0.14 


0.14 


0.14 


0.14 


0.14 


0.14 


0.14 


0.14 


2 


2-^ 


0.16 


0.16 


0.16 


0.16 


0.16 


0.15 


0.15 


0.16 


0.16 


0.16 


0.16 


3 


2-4 


0.16 


0.16 


0.16 


0.16 


0.16 


0.16 


0.16 


0.16 


0.16 


0.16 


0.16 



Table 7.2. Norm of the error propagator E 
to IIPG discretization with a = 4a*. 



{L — Ag^ A) for A corresponding 



More experiments for the IP(/3)-l methods can be found in [9]. 
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Appendix A. Construction of an Interpolation Operator 
We now construct an interpolation operator which satisfies the approximation and stability 



properties (5.15)-(5.16) in Lemma 5.5 To begin with, let us introduce some notation. Given a 
conforming triangulation Th, recall that £h is the set of edges/faces of Th- Let Sh C H^{yi) be 
the conforming P'^(7h) nC°(J]) La grange finite element space (quadratics in d = 2 and cubics for 
d = 3). We split the set of DOFs of Sh into two subsets C{Th) and M{Th), where C{Th) contains 
the DOFs of corresponding to the barycenters of the edges/faces in S^, and M{Th) contains 
all the remaining DOFs of Sh. We also denote the restriction of Th, £hi ^h^ or Sh on a given 
subdomain G by Th{G), £h{G), V^'^{G) or Sh{G), respectively. 

Let 7^ be a coarser mesh, i.e., Th is either the same as 7^ or a refinement of it with h < h. We 
now start building the operator : — t- ^°"^) where we recall that is the piecewise 



conforming finite element space defined on 7^. The basic idea is to, in each subdomain Oj, embed 

/cc 
h 

operator. 



V^^{^li) into Shi^i). Then we interpolate the result in V~°'^^ on 7^ using a quasi-interpolation 
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To embed into Sh, we modify the inclusion operator introduced in [16], and define it at 
the subdomain level as follows. For any v G we define Ei : V^'^{^i) — )• Sh{^i) on each 

subdomain $7j as: 

' v{p), ifp£C{Th)nUi 
#w; T.T&M^ ^t{p), if p g (Tft) n Qi (a.i) 



{Eiv){p) 



#M^EeeAfa^e(p)' H p £ ^f (Jh) H d^i, 



where Mp := {T G Th{^i) ■ p G dT} is the set of elements sharing p and Mp := {e G £h{^i) '■ 
e C d^i s.t. p G de} is the set of edges on dQi sharing the DOF p. Here #Mp and #Mp 
are the cardinality of these sets respectively, and vt, Ve are the restriction of u on T and e 
respectively. 

Observe that, this construction differs from the one in [161 Equation (3.1)] in the treatment of 



the DOFs on dfli. From (A.I), for each DOF p G d^i, {Eiv){p) contains only the contributions 
of V from the boundary of Vti, not from the interior. Therefore, it is obvious that {Eiv){p) = 
{Ejv){p) for any DOF p in the interior of the interface F = dVti n d^j {i ^ j) between the 
subdomains fij and Vtj. This special treatment at boundary points guarantees that the global 
function r^l^. := EiV is continuous for the points in the interior of each interface. However, this 
global function r] will generally be multi- valued at the points on dV. 



Although the construction of Ei in (A.I) is different from [16], the same analysis in |T6j can 



be carried out here. We summarize the the properties of Ei below, and omit the detailed proof. 

Lemma A.I. The linear operator Ei : V^'^{^li) — )• Sh{^i) defined in (A.I) satisfies that 

\EMi,n,, - bli./iA, and \\v - EivW^^^^ < h \v\^^h^^^ , G T/f ^. 

Let Qi : H^{^i) Vf°'^^{ni) and Qr : ^^F) ^ Vf^^'^T) be the Scott-Zhang quasi- 
interpolation operators on J7j and on the interface F C ilj, respectively. We now recall the 
definition and main properties of these operators. In the sequel, we should denote a generic 
vertex of 7^ by p. Let ojp := IJ ^ 'Tjii^i) '■ P ^ C Qi be the local patch containing 
p, and LOT '■= U {'^p • P ^ 9T} for each T G Tj^{Qi). Similarly, on the interface F, we define 
Op := U{e G : e C F and G ae} C F and Oe := [jWp ■ P & de} for each e G £-f^{T). For any 
vertex p, let 4>p G V-""^^ be the nodal basis function, and define the dual basis 9p G V~°'^^{ujp) 
such that 



9pvdx = v{p), Vw G Vf^. 



To define Qi, let us choose somen T C ujp for each vertex p G Tj^{Qi). Then, the Scott-Zhang 
quasi-interpolation operator is denned by 

QiV= Yl (fdpVdxj^p, rjeH\Qi). 

The operator Qr is defined similarly, but restricted on the interface F. Both operators enjoy 
the following approximation and stability properties (see }49| [53] for a proof): 

Lemma A. 2. For any rj G H^{VLi), the operator Qi : H^{Qi) — )• V~°^^{Qi) satisfies: 

||Q^??||o,T< ||r?||o,^^, \\Q^v\\l,T < , II " Qi)^llo,T < /^hHi.c.^, VT G r^(17i). (A.2) 

For any ^ G H^{r), the operator Qr : H^{r) — )• V~°°^(F) satisfies the following properties: 

IIQr?llo,e<||^||o,Oe, IIU-Qr)e)l|o,e</i||elli,o. VeG£:^(F). (A.3) 

Note that the choice of T may not be unique. 
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Furthermore, both operators are linear preserving; i.e. QiT] = rj for any rj G V~"'^^{^i), anc 
similarly QrC = ^ for any £^ G V~°^^(T). 

Now we are ready to define the interpolation operator Pj^ : — )• V~°'^^: 



{QiE,v){p), ifpea 
Ph"^ ]\nM = { {QvEiv) (p), if p G int(r) for eacli side T C dQi , 
0, elsewliere 



(A.4) 



where int(r) is the interior of T. From the definition of Ei in ( |A.l ), if p is a vertex of 7^ in the 
interior of the interface T = iliCifij, we have {Eiv){p) = {Ejv){p), which implies (QrEiv) (p) = 
(QrEjv) (p). The special treatment for the interface in ( |A.4 ) guarantees the global continuity 

of Pj^v. Thus, Pj^v G ^""^^ is well-defined. Now, we show that the operator defined in (A.4) 
does satisfy the approximation and stability properties (5.15)-(5.16): 

Lemma A. 3. For any v G V^^, the operator : — )• V~°^^ satisfies 



\\{I-Phv\\o,. < h\logi2h/h)\'/^vh,H,n, 



(A.5) 
(A.6) 



Proof. The proof follows the ideas from [14', Lemma 4.6], adapted to the present situation. We 
start by showing (A.5). Using triangle inequality. Lemma A.l, together with the approximation 
result (A.2) of the Qi from Lemma A. 2, we have 

\\v - Ptv\Wn, < h - QiEivWo^n, + \\QiEiV - pMo^^^ 

< \\v - Eiv\\o,n, + - Qi)Eiv\\o,n, + WQiEtv - P^v\\o,nt 

< h\v\i^h,n, +h\\Eiv\\i^n^ + \\QiEiV - P^v||o,n, 



< ^bli./iA + hMi^h,n, + WQiEiV - Pj^vl 



(A.7) 



Hence, to show the inequality (A.5) we only need to estimate \\QiEiV — P/*f ||o,f2i- 

To simplify the notation, throughout the proof we set X = P~v G as defined in (|X4|, 

and denote Xi •= QiEiV. From the definition of P^ in (A.4), x{p) = Xi{p) when p is a vertex 
of 7^ in the interior of fli, and they are different only on the boundary vertices. So by using 
discrete norm, we have 



\QiE^v-Pi'v 



hV\\o,n, 



rcen, per 



Vcdni \pGint(r) pedv 



^ E E h\\QrE,v-x^\\le + h''\\x 



illo,9r 



Below, we try to bound those two terms appearing in the last expression of (A.8). 
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For the first term in (A. 8), we observe that Xi = QrXi by Lemma A. 2, Then by the L 
stabihty property (A. 3) of Qr, we obtain 
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\QrE^v-x^ 



i\\0,e 



hWQrEiV - QrXiWle < h\\EiV - x^\ 



0,Oe 



< hih)-'\\{I - Qi)Eiv\\l^^ + h'\EiV - QiEMU^ 

where in the second inequahty, we used the standard trace inequahty (cf. [HI Lemma 2.1]), and 
in the last step we used the properties (A. 2) of Qj. Here := U{T E 7^(^^j) : dT r\ Oe 7^ 0}. 
Summing up the above inequality for all edges/faces on we obtain that 

To bound the second term in (A. 8) we have to distinguish between the 2D and 2>D cases. In 
the 2D case, P is a one-dimensional edge of rij, so dV reduces to its two endpoints, say {p, (?}. 
Hence, 

llXillo.ar = {\Xi{p)? + \Xi{(l)?) < IIXi|lo,oo,.^p + llXillo.oo.u;, , dT = {p,q} . 

To bound each of the above two terms on the right side, we use the two-dimensional discrete 
Sobolev inequality Lemma 2.3]; 



IXi||0,( 



< C ( log ^^""^^'"-^ 



h 



1/2 



\Xi\\ i-,ujp 



(A.IO) 



So summing over all F C d^li the resulting estimate, we finally get 



E 

rcdUi 



rcdfii pear 
'2h\ , 



/ diam(a;p 



IXiWlcOr. ~ log 



2h 
~h 



I All 1 1,0,: 



log 



h 



QiE^v\\ln, < log 



2h 
~h 



(A.ll) 



where in the second inequality we used the fact diam(a;p) ~ 2h, and in the last step we used 



the inequality (A.2) of Qi and the properties of Ei in Lemma A.l 



In 3D, F C dQi is a two-dimensional face of So 9F is a union of edges in the triangulation 
{e € Sj^ : e C dT}. In this case, we use the following the discrete Sobolev inequality [141 Lemma 
2.4] (instead of (|A.10[) in 2D case): 



\Xi\\o,dr 



IIXi|lo,e ^ 



eCdr 



/ diam(a;e 



lvl|2 



Summing the above estimate over all F C 90j and using, as before, the inequality (A.2) of Qi 
together with the properties of Ei given in Lemma |A.1[ we find 

/ diam(a;e 



rcdcii 



'2h\ , 



\Xi\\l,u, < log 



2h 



E El 

rc9Qi ecar 



vl|2 



^ log 



h 



XilliA ~ l°g 



2h 
~h 



\l,h,Qi 



(A.12) 



Now, substituting (A.12) (or (A.ll) when d = 2) and (A.9) into (A.8), we finally get 



\QiEiV - PlivWl^^ = \\x - XilloA < f^^\Hlh,n, + log 



2h 



\l,h,Qi 



30 



B. AYUSO DE DIGS, M. HOLST, Y. ZHU, AND L. ZIKATANOV 



The inequality (A. 5) then follows by inserting the above estimate in (A. 7). 

Finally we show the stability of p} (jXB]). Note that pj^v £ yi^°^^and v £ V^^. To deal 
with possibly different mesh sizes we consider the local L^-projection Vt ■ L'^iT) — > F^{T) for 
any T £ Ij^. For h > h, such an element is the union of other subelements in the partition Th- 
Then, adding and subtracting Vtv, triangle inequality together with inverse inequality and the 



approximation property (A.5), gives 



\Phv\i,T < \PhV - Vtv\i,t + \Ptv\i,t < C{hr'^\\Pliv -'Ptv\\o,t + \Ptv\i,t 
< C(h)-^ (wpI^v - v\\o,T + \\v - Vtv\\o,t) + C\v\i,T 



< C{h)-^\\Pliv - vWo^T + C\\v\\i,T ■ 

The Stability now follows immediately, by summing over all elements T G i^i, using the definition 
of the weighted ff^-semi-norm and the weighted L^-norm together with the approximation result 
already shown: 



\Ptv 



< Ch 



-1 



■zn ^''^ 

< Ch'^h ( log ( — 1 1 + Iblli.^K.n 



< 



log 



2h 




and the proof is complete. 



□ 
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